Based on the alpha and beta diversity results, I am not expecting to find a huge difference between groups in terms of bacterial composition. Please note that many of these code are from the Microbiome Intensive Course (MIC Course, Duke University, 2023) and thus the work of Dr. Pixu Shi.
This is after talking with Dr. Shi about what analyses I want to run.
First, update the ANCOMBC package if necessary.
#if (!requireNamespace("BiocManager", quietly=TRUE))
# install.packages("BiocManager")
#BiocManager::install("ANCOMBC")
# rm(list = ls(all = TRUE)) # Unload all packages if necessary
library(phyloseq)
library(fs)
library(rstatix)
library(ggpubr)
library(microViz)
library(ANCOMBC)
library(pals)
library(RColorBrewer)
library(colorBlindness)
library(microbiome)
library(ggrepel)
library(DT)
library(plotly)
library(DESeq2)
library(parallel)
library(dplyr)
library(tidyr)
library(readr)
library(tibble)
library(stringr)
library(cowplot)
rm(list=ls()) # clear environment
git.dir = file.path("/hpc/group/dmc/Eva/GLP_KO_Microbiome_final") # CHANGE ME
ps.rds <- file.path(git.dir, "ps.GLP1.rds")
ps <- read_rds(ps.rds)
ps
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
fig.dir = file.path(git.dir, "figures")
map.file = file.path(git.dir, "Metadata_GLP1.csv")
meta.df = read_csv(map.file, show_col_types = FALSE)
I will now aggregate at the genus level, as it is the smallest taxonomic rank that analyses will be run on. In other words, I do not plan on running differential abundance analysis at ASV level.
Let’s prepare the phyloseq object for differential abundance analysis (DAA). Since we are interested in 1) differentially abundant taxa and 2) correlation analysis, we will limit the taxa to only ones that appear in > 3 samples for each sample type.
To avoid running analyses on extremely rare taxa that may be sequencing artifacts, the threshold for minimum number of counts is set at 10. This is, of course, arbitrary. 10 was chosen as 1% of number of read threshold of 1000 above for pruning samples. Based on how I wrote the code, taxa must be present have > 10 reads or more in > 3 samples.
Furthermore, categorical variables should be re-leveled as factors.
# Subset to true fecal samples only and remove tree
ps %>% subset_samples(Sample_type == "Fecal" & Type == "True Sample") -> ps_fecal
ps_fecal_notree <- phyloseq(otu_table(ps_fecal), tax_table(ps_fecal), sample_data(ps_fecal))
# Categorical variables should be factors
ps_fecal_notree@sam_data$Genotype <- factor(ps_fecal_notree@sam_data$Genotype, levels = c("WT", "KO"))
ps_fecal_notree@sam_data$Sex <- factor(ps_fecal_notree@sam_data$Sex, levels = c("Female", "Male"))
ps_fecal_notree@sam_data$Intranasal_Treatment <- factor(ps_fecal_notree@sam_data$Intranasal_Treatment, levels = c("PBS", "HDM"))
ps_fecal_notree@sam_data$Surgery <- factor(ps_fecal_notree@sam_data$Surgery, levels = c("None", "Sham", "VSG"))
ps_fecal_notree@sam_data$Mouse %>% unique() -> mouse_levels
ps_fecal_notree@sam_data$Mouse <- factor(ps_fecal_notree@sam_data$Mouse, levels = mouse_levels)
ps_fecal_notree@sam_data$Diet <- factor(ps_fecal_notree@sam_data$Diet, levels = c("Normal_Chow", "High_Fat_Diet"))
ps_fecal_notree@sam_data$Group <- factor(ps_fecal_notree@sam_data$Group, levels = c("NA", "Control", "1", "2"))
ps_fecal_notree@sam_data$Week <- factor(ps_fecal_notree@sam_data$Week, levels = c("0", "10", "13"))
ps_fecal_notree
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
# Aggregate at genus level
ps_g_notree_fecal <- tax_glom(ps_fecal_notree, "Genus")
taxa_names(ps_g_notree_fecal) <- tax_table(ps_g_notree_fecal)[, 'Genus']
ps_g_notree_fecal
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 314 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 314 taxa by 6 taxonomic ranks ]
We will prune the phyloseq object at each phylogenetic rank at and above genus. Since we plan on running correlation analysis after differential abundance, let’s test for taxa that are present in at least 20% of all samples (arbitrary) with greater than 10 reads in each sample (also arbitrary. It’s also 1% of 1000 reads, which was the threshold used to remove samples from analysis).
x` Genus
# Pruning: set threshold number
threshold_num = 0.20 * nrow(sample_data(ps_fecal))
nreads_prune = 10
# Prune taxa at each phylogenetic rank
ps_g_notree_fecal %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_genus.prune
ps_genus.prune
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 56 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 56 taxa by 6 taxonomic ranks ]
ps_g_notree_fecal %>% transform_sample_counts(function(x) x/sum(x)) -> ps_genus.ts
ps_genus.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 314 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 314 taxa by 6 taxonomic ranks ]
subset_taxa(ps_genus.ts, ps_genus.ts@tax_table[, "Genus"] %in% ps_genus.prune@tax_table[, "Genus"]) -> ps_genus.prune.ts
ps_genus.prune.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 56 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 56 taxa by 6 taxonomic ranks ]
rowSums(ps_genus.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> genus_prune_RelAb
summary(genus_prune_RelAb$RelAb)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 67.11 99.32 99.65 99.25 99.85 100.00
ps_fecal_notree %>% tax_glom("Family") -> ps_Family
taxa_names(ps_Family) <- tax_table(ps_Family)[, 'Family']
ps_Family
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 169 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 169 taxa by 6 taxonomic ranks ]
ps_Family %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Family.ts
ps_Family.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 169 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 169 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank
ps_Family %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Family.prune
ps_Family.prune
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 28 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 28 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Family.ts, ps_Family.ts@tax_table[, "Family"] %in% ps_Family.prune@tax_table[, "Family"]) -> ps_Family.prune.ts
ps_Family.prune.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 28 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 28 taxa by 6 taxonomic ranks ]
rowSums(ps_Family.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Family_prune_RelAb
summary(Family_prune_RelAb)
## RelAb
## Min. : 71.60
## 1st Qu.: 99.68
## Median : 99.85
## Mean : 99.54
## 3rd Qu.: 99.94
## Max. :100.00
ps_fecal_notree %>% tax_glom("Order") -> ps_Order
taxa_names(ps_Order) <- tax_table(ps_Order)[, 'Order']
ps_Order
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 109 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 109 taxa by 6 taxonomic ranks ]
ps_Order %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Order.ts
ps_Order.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 109 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 109 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank
ps_Order %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Order.prune
ps_Order.prune
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 18 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 18 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Order.ts, ps_Order.ts@tax_table[, "Order"] %in% ps_Order.prune@tax_table[, "Order"]) -> ps_Order.prune.ts
ps_Order.prune.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 18 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 18 taxa by 6 taxonomic ranks ]
rowSums(ps_Order.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Order_prune_RelAb
summary(Order_prune_RelAb)
## RelAb
## Min. : 71.63
## 1st Qu.: 99.76
## Median : 99.89
## Mean : 99.56
## 3rd Qu.: 99.96
## Max. :100.00
ps_fecal_notree %>% tax_glom("Class") -> ps_Class
taxa_names(ps_Class) <- tax_table(ps_Class)[, 'Class']
ps_Class
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 42 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 42 taxa by 6 taxonomic ranks ]
ps_Class %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Class.ts
ps_Class.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 42 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 42 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank
ps_Class %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Class.prune
ps_Class.prune
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 12 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 12 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Class.ts, ps_Class.ts@tax_table[, "Class"] %in% ps_Class.prune@tax_table[, "Class"]) -> ps_Class.prune.ts
ps_Class.prune.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 12 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 12 taxa by 6 taxonomic ranks ]
rowSums(ps_Class.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Class_prune_RelAb
summary(Class_prune_RelAb)
## RelAb
## Min. : 99.99
## 1st Qu.:100.00
## Median :100.00
## Mean :100.00
## 3rd Qu.:100.00
## Max. :100.00
ps_fecal_notree %>% tax_glom("Phylum") -> ps_Phylum
taxa_names(ps_Phylum) <- tax_table(ps_Phylum)[, 'Phylum']
ps_Phylum
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
ps_Phylum %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Phylum.ts
ps_Phylum.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank
ps_Phylum %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Phylum.prune
ps_Phylum.prune
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Phylum.ts, ps_Phylum.ts@tax_table[, "Phylum"] %in% ps_Phylum.prune@tax_table[, "Phylum"]) -> ps_Phylum.prune.ts
ps_Phylum.prune.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 6 taxonomic ranks ]
rowSums(ps_Phylum.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Phylum_prune_RelAb
summary(Phylum_prune_RelAb)
## RelAb
## Min. :100
## 1st Qu.:100
## Median :100
## Mean :100
## 3rd Qu.:100
## Max. :100
Let’s visualize at Phylum level to see how much information was kept:
plot_bar(ps_genus.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Genus", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps)
plot_bar(ps_Family.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Family", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps)
plot_bar(ps_Order.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Order", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps)
plot_bar(ps_Class.prune.ts, x = "Label", fill="Class") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Class, fill=Class), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Class", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps)
plot_bar(ps_Phylum.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Phylum", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps)
meta.df %>% filter(Type == "True Sample" & Sample_type == "Fecal") %>%
filter(Label %in% phyloseq::sample_data(ps_g_notree_fecal)$Label) -> meta.fecal
meta.fecal
Let’s see how many mice we have:
meta.fecal %>% filter(Diet == "High_Fat_Diet" & Week == "13") -> meta.f.hfd
meta.f.hfd %>%
dplyr::count(Surgery, Genotype) %>%
arrange(n)
meta.f.hfd %>%
dplyr::count(Surgery) %>%
arrange(n)
nrow(meta.f.hfd)
## [1] 45
meta.f.hfd %>% distinct(Mouse) %>% nrow()
## [1] 45
ps_genus.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Genus
tse.fecal.Q1.Genus
## class: TreeSummarizedExperiment
## dim: 56 45
## metadata(0):
## assays(1): counts
## rownames(56): Clostridium sensu stricto 1 Anaerotruncus ... Romboutsia
## Bifidobacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Family.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Family
tse.fecal.Q1.Family
## class: TreeSummarizedExperiment
## dim: 28 45
## metadata(0):
## assays(1): counts
## rownames(28): Peptococcaceae Clostridiaceae ... Anaerovoracaceae
## Bifidobacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Order.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Order
tse.fecal.Q1.Order
## class: TreeSummarizedExperiment
## dim: 18 45
## metadata(0):
## assays(1): counts
## rownames(18): Peptococcales Clostridiales ...
## Peptostreptococcales-Tissierellales Bifidobacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Class.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Class
tse.fecal.Q1.Class
## class: TreeSummarizedExperiment
## dim: 12 45
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
## Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Phylum.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Phylum
tse.fecal.Q1.Phylum
## class: TreeSummarizedExperiment
## dim: 10 45
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
## Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
set.seed(123) # Set seed for reproducibility
n_cl_requested = 24
# Genus level
fecal.Q1.Genus <- ANCOMBC::ancom(data=tse.fecal.Q1.Genus, assay_name="counts", tax_level="Genus",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Genus = fecal.Q1.Genus$res
res.fecal.Q1.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Genus.sig
tab.fecal.Q1.Genus.sig
set.seed(123) # Set seed for reproducibility
# Family level
fecal.Q1.Family <- ANCOMBC::ancom(data=tse.fecal.Q1.Family, assay_name="counts", tax_level="Family",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Family = fecal.Q1.Family$res
res.fecal.Q1.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Family.sig
tab.fecal.Q1.Family.sig
set.seed(123) # Set seed for reproducibility
# Order level
fecal.Q1.Order <- ANCOMBC::ancom(data=tse.fecal.Q1.Order, assay_name="counts", tax_level="Order",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Order = fecal.Q1.Order$res
res.fecal.Q1.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Order.sig
tab.fecal.Q1.Order.sig
set.seed(123) # Set seed for reproducibility
# Class level
fecal.Q1.Class <- ANCOMBC::ancom(data=tse.fecal.Q1.Class, assay_name="counts", tax_level="Class",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Class = fecal.Q1.Class$res
res.fecal.Q1.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Class.sig
tab.fecal.Q1.Class.sig
set.seed(123) # Set seed for reproducibility
# Phylum level
fecal.Q1.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q1.Phylum, assay_name="counts", tax_level="Phylum",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Phylum = fecal.Q1.Phylum$res
res.fecal.Q1.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Phylum.sig
tab.fecal.Q1.Phylum.sig
tab.fecal.Q1.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q1.Family.sig$Tax_rank <- "Family"
tab.fecal.Q1.Order.sig$Tax_rank <- "Order"
tab.fecal.Q1.Class.sig$Tax_rank <- "Class"
tab.fecal.Q1.Phylum.sig$Tax_rank <- "Phylum"
rbind(tab.fecal.Q1.Genus.sig,
tab.fecal.Q1.Family.sig,
tab.fecal.Q1.Order.sig,
tab.fecal.Q1.Class.sig,
tab.fecal.Q1.Phylum.sig) -> tab.fecal.Q1
tab.fecal.Q1$Test <- "ANCOM II, fecal microbiome: all HFD, wk 10 & 13. Structural zeroes from Surgery, adjust with Genotype"
tab.fecal.Q1
ds.fecal.Q1_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q1.Genus, ~ Surgery * Genotype)
## converting counts to integer mode
dds.fecal.Q1_Genus <- DESeq2::DESeq(ds.fecal.Q1_Genus)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Genus)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 580367 2515 33645 12845 12897. 5729.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 56 0 37.4
resultsNames(dds.fecal.Q1_Genus)
## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Genotype_KO_vs_WT" "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"
results(dds.fecal.Q1_Genus, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
filter(padj < 0.05)
results(dds.fecal.Q1_Genus, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q1.Genus.Res
fecal.Q1.Genus.Res$Comparison <- c("Surgery: VSG over None")
results(dds.fecal.Q1_Genus, name = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Genus, name = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Sham surgery * KO genotype")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res
results(dds.fecal.Q1_Genus, name = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("VSG * KO genotype")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res
fecal.Q1.Genus.Res$Compare_taxon <- "Genus"
fecal.Q1.Genus.Res
fecal.Q1.Genus.Res %>% filter(row %in% tab.fecal.Q1.Genus.sig)
ds.fecal.Q1_Family <- DESeq2::DESeqDataSet(tse.fecal.Q1.Family, ~ Surgery * Genotype)
## converting counts to integer mode
dds.fecal.Q1_Family <- DESeq2::DESeq(ds.fecal.Q1_Family)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 2 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Family)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 721636 3576 43678 15928 16036. 7515.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 28 0 21.3
resultsNames(dds.fecal.Q1_Family)
## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Genotype_KO_vs_WT" "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"
results(dds.fecal.Q1_Family, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Family.Res
fecal.Q1.Family.Res$Comparison <- c("Surgery: Sham over None")
results(dds.fecal.Q1_Family, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Family, name = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Family, name = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Family, name = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05)
fecal.Q1.Family.Res$Compare_taxon <- "Family"
fecal.Q1.Family.Res
fecal.Q1.Family.Res %>% filter(row %in% tab.fecal.Q1.Family.sig$taxon)
ds.fecal.Q1_Order <- DESeq2::DESeqDataSet(tse.fecal.Q1.Order, ~ Surgery * Genotype)
## converting counts to integer mode
dds.fecal.Q1_Order <- DESeq2::DESeq(ds.fecal.Q1_Order)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Order)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 724646 3587 43755 15943 16103. 7577.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 18 0 13.9
resultsNames(dds.fecal.Q1_Order)
## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Genotype_KO_vs_WT" "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"
results(dds.fecal.Q1_Order, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Order.Res
fecal.Q1.Order.Res$Comparison <- c("Surgery:S Sham over none")
results(dds.fecal.Q1_Order, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
results(dds.fecal.Q1_Order, name = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Order, name = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Order, name = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05)
fecal.Q1.Order.Res$Compare_taxon <- "Order"
fecal.Q1.Order.Res
fecal.Q1.Order.Res %>% filter(row %in% tab.fecal.Q1.Order.sig$taxon)
ds.fecal.Q1_Class <- DESeq2::DESeqDataSet(tse.fecal.Q1.Class, ~ Surgery * Genotype)
## converting counts to integer mode
dds.fecal.Q1_Class <- DESeq2::DESeq(ds.fecal.Q1_Class)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 2 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Class)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 729737 3587 43771 15951 16216. 7521.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 12 0 9.13
resultsNames(dds.fecal.Q1_Class)
## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Genotype_KO_vs_WT" "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"
results(dds.fecal.Q1_Class, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
filter(padj < 0.05)
results(dds.fecal.Q1_Class, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q1.Class.Res
fecal.Q1.Class.Res$Comparison <- c("Surgery: VSG over None")
results(dds.fecal.Q1_Class, name = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Class, name = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Class, name = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05)
fecal.Q1.Class.Res$Compare_taxon <- "Class"
fecal.Q1.Class.Res
tab.fecal.Q1.Class.sig
ds.fecal.Q1_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q1.Phylum, ~ Surgery * Genotype)
## converting counts to integer mode
dds.fecal.Q1_Phylum <- DESeq2::DESeq(ds.fecal.Q1_Phylum)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Phylum)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 729772 3587 43771 15953 16217. 7521.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 10 0 7.84
resultsNames(dds.fecal.Q1_Phylum)
## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Genotype_KO_vs_WT" "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"
results(dds.fecal.Q1_Phylum, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Phylum.Res
fecal.Q1.Phylum.Res$Comparison <- c("Surgery: Sham over None")
results(dds.fecal.Q1_Phylum, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res
results(dds.fecal.Q1_Phylum, name = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q1_Phylum, name = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Sham surgery * KO genotype")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res
results(dds.fecal.Q1_Phylum, name = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
filter(padj < 0.05)
fecal.Q1.Phylum.Res$Compare_taxon <- "Phylum"
fecal.Q1.Phylum.Res
tab.fecal.Q1.Phylum.sig
rbind(fecal.Q1.Genus.Res,
fecal.Q1.Family.Res,
fecal.Q1.Class.Res,
fecal.Q1.Order.Res,
fecal.Q1.Phylum.Res) -> fecal.Q1.deseq
fecal.Q1.deseq$Test <- "DESeq2, Fecal microbiome: all HFD, wk 10 + 13, Surgery * Genotype"
fecal.Q1.deseq
tab.fecal.Q1 %>% filter(taxon %in% fecal.Q1.deseq$row) -> ancom.intersect.hfd
ancom.intersect.hfd
fecal.Q1.deseq %>% filter(row %in% tab.fecal.Q1$taxon)
psmelt(ps_genus.prune.ts) -> Genus.prune.melt
psmelt(ps_Family.prune.ts) -> Family.prune.melt
psmelt(ps_Order.prune.ts) -> Order.prune.melt
psmelt(ps_Class.prune.ts) -> Class.prune.melt
psmelt(ps_Phylum.prune.ts) -> Phylum.prune.melt
bind_rows(
Genus.prune.melt,
Family.prune.melt,
Order.prune.melt,
Class.prune.melt,
Phylum.prune.melt
) -> fecal.RA.Full
fecal.RA.Full$Surgery <- factor(fecal.RA.Full$Surgery, levels = c("None", "Sham", "VSG"))
fecal.RA.Full %>%
dplyr::select(OTU:Label, Sample_type:Mouse, Genotype:Surgery, Week) -> fecal.RA.Full
head(fecal.RA.Full)
tail(fecal.RA.Full)
fecal.RA.Full %>% filter(Diet == "High_Fat_Diet" & Week != "0") -> fecal.graph.hfd
head(fecal.graph.hfd)
tail(fecal.graph.hfd)
Let’s see how many mice we have:
meta.fecal %>% filter(Surgery == "None") -> meta.f.nosurgery
meta.f.nosurgery %>%
dplyr::count(Week, Diet, Genotype) %>%
arrange(n)
meta.f.nosurgery %>%
dplyr::count(Week, Diet) %>%
arrange(n)
I will run Diet * Week + Genotype.
meta.f.nosurgery %>% distinct(Mouse) %>% nrow()
## [1] 58
ps_genus.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Genus
tse.fecal.Q2.Genus
## class: TreeSummarizedExperiment
## dim: 56 120
## metadata(0):
## assays(1): counts
## rownames(56): Clostridium sensu stricto 1 Anaerotruncus ... Romboutsia
## Bifidobacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Family.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Family
tse.fecal.Q2.Family
## class: TreeSummarizedExperiment
## dim: 28 120
## metadata(0):
## assays(1): counts
## rownames(28): Peptococcaceae Clostridiaceae ... Anaerovoracaceae
## Bifidobacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Order.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Order
tse.fecal.Q2.Order
## class: TreeSummarizedExperiment
## dim: 18 120
## metadata(0):
## assays(1): counts
## rownames(18): Peptococcales Clostridiales ...
## Peptostreptococcales-Tissierellales Bifidobacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Class.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Class
tse.fecal.Q2.Class
## class: TreeSummarizedExperiment
## dim: 12 120
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
## Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Phylum.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Phylum
tse.fecal.Q2.Phylum
## class: TreeSummarizedExperiment
## dim: 10 120
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
## Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
set.seed(123) # Set seed for reproducibility
# Genus level
fecal.Q2.Genus <- ANCOMBC::ancom(data=tse.fecal.Q2.Genus, assay_name="counts", tax_level="Genus",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Genus = fecal.Q2.Genus$res
res.fecal.Q2.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Genus.sig
tab.fecal.Q2.Genus.sig
set.seed(123) # Set seed for reproducibility
# Family level
fecal.Q2.Family <- ANCOMBC::ancom(data=tse.fecal.Q2.Family, assay_name="counts", tax_level="Family",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Family = fecal.Q2.Family$res
res.fecal.Q2.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Family.sig
tab.fecal.Q2.Family.sig
set.seed(123) # Set seed for reproducibility
# Order level
fecal.Q2.Order <- ANCOMBC::ancom(data=tse.fecal.Q2.Order, assay_name="counts", tax_level="Order",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Order = fecal.Q2.Order$res
res.fecal.Q2.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Order.sig
tab.fecal.Q2.Order.sig
set.seed(123) # Set seed for reproducibility
# Class level
fecal.Q2.Class <- ANCOMBC::ancom(data=tse.fecal.Q2.Class, assay_name="counts", tax_level="Class",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Class = fecal.Q2.Class$res
res.fecal.Q2.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Class.sig
tab.fecal.Q2.Class.sig
set.seed(123) # Set seed for reproducibility
# Phylum level
fecal.Q2.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q2.Phylum, assay_name="counts", tax_level="Phylum",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested)
## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Phylum = fecal.Q2.Phylum$res
res.fecal.Q2.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Phylum.sig
tab.fecal.Q2.Phylum.sig
tab.fecal.Q2.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q2.Family.sig$Tax_rank <- "Family"
tab.fecal.Q2.Order.sig$Tax_rank <- "Order"
tab.fecal.Q2.Class.sig$Tax_rank <- "Class"
tab.fecal.Q2.Phylum.sig$Tax_rank <- "Phylum"
rbind(tab.fecal.Q2.Genus.sig,
tab.fecal.Q2.Family.sig,
tab.fecal.Q2.Order.sig,
tab.fecal.Q2.Class.sig,
tab.fecal.Q2.Phylum.sig) -> tab.fecal.Q2
tab.fecal.Q2$Test <- "ANCOM II, fecal microbiome: no surgery, main effect Diet (structural zeroes) and week + genotype"
tab.fecal.Q2
ds.fecal.Q2_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q2.Genus, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Genus <- DESeq2::DESeq(ds.fecal.Q2_Genus)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Genus)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1231191 2007 45648 9976. 10260. 5536.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 56 0 40.0
resultsNames(dds.fecal.Q2_Genus)
## [1] "Intercept" "Diet_High_Fat_Diet_vs_Normal_Chow"
## [3] "Week_10_vs_0" "Week_13_vs_0"
## [5] "Genotype_KO_vs_WT" "DietHigh_Fat_Diet.Week10"
## [7] "DietHigh_Fat_Diet.Week13"
results(dds.fecal.Q2_Genus, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05)
results(dds.fecal.Q2_Genus, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Genus.Res
fecal.Q2.Genus.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
results(dds.fecal.Q2_Genus, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
results(dds.fecal.Q2_Genus, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
fecal.Q2.Genus.Res$Compare_taxon <- "Genus"
fecal.Q2.Genus.Res
ds.fecal.Q2_Family <- DESeq2::DESeqDataSet(tse.fecal.Q2.Family, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Family <- DESeq2::DESeq(ds.fecal.Q2_Family)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Family)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1566638 2590 55908 12286. 13055. 6960.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 28 0 23.2
results(dds.fecal.Q2_Family, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05)
results(dds.fecal.Q2_Family, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Family.Res
fecal.Q2.Family.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
results(dds.fecal.Q2_Family, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
results(dds.fecal.Q2_Family, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
fecal.Q2.Family.Res$Compare_taxon <- "Family"
fecal.Q2.Family.Res
ds.fecal.Q2_Order <- DESeq2::DESeqDataSet(tse.fecal.Q2.Order, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Order <- DESeq2::DESeq(ds.fecal.Q2_Order)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Order)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1589033 2592 56999 12546. 13242. 7087.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 18 0 16.0
results(dds.fecal.Q2_Order, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05)
results(dds.fecal.Q2_Order, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Order.Res
fecal.Q2.Order.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res
results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res
results(dds.fecal.Q2_Order, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q2_Order, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res
fecal.Q2.Order.Res$Compare_taxon <- "Order"
fecal.Q2.Order.Res
ds.fecal.Q2_Class <- DESeq2::DESeqDataSet(tse.fecal.Q2.Class, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Class <- DESeq2::DESeq(ds.fecal.Q2_Class)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Class)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1594382 2595 57079 12562 13287. 7112.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 12 0 10.6
results(dds.fecal.Q2_Class, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05)
results(dds.fecal.Q2_Class, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Class.Res
fecal.Q2.Class.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res
results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res
results(dds.fecal.Q2_Class, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q2_Class, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res
fecal.Q2.Class.Res$Compare_taxon <- "Class"
fecal.Q2.Class.Res
ds.fecal.Q2_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q2.Phylum, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Phylum <- DESeq2::DESeq(ds.fecal.Q2_Phylum)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Phylum)
## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1594436 2595 57083 12562 13287. 7112.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 10 0 8.79
results(dds.fecal.Q2_Phylum, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05)
results(dds.fecal.Q2_Phylum, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Phylum.Res
fecal.Q2.Phylum.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res
results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res
results(dds.fecal.Q2_Phylum, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05)
results(dds.fecal.Q2_Phylum, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res
fecal.Q2.Phylum.Res$Compare_taxon <- "Phylum"
fecal.Q2.Phylum.Res
rbind(fecal.Q2.Genus.Res,
fecal.Q2.Family.Res,
fecal.Q2.Class.Res,
fecal.Q2.Order.Res,
fecal.Q2.Phylum.Res) -> fecal.Q2.deseq
fecal.Q2.deseq$Test <- "DESeq2, Fecal microbiome: no surgery, ~ Diet * Week + Genotype"
fecal.Q2.deseq
tab.fecal.Q2 %>% filter(taxon %in% fecal.Q2.deseq$row) -> ancom.intersect.nosurgery
ancom.intersect.nosurgery
ancom.intersect.nosurgery %>% arrange(Tax_rank, taxon)
fecal.Q2.deseq %>% filter(row %in% tab.fecal.Q2$taxon) %>%
dplyr::select(row, Comparison, Compare_taxon, Test) %>%
group_by(row) %>%
mutate(n.ind = paste0(1:n())) %>%
pivot_wider(
names_from = "n.ind",
values_from = Comparison
) -> fecal.Q2.deseq.intersect.wide
fecal.Q2.deseq.intersect.wide
fecal.RA.Full %>% filter(Surgery == "None") -> fecal.graph.nosurgery
head(fecal.graph.nosurgery)
tail(fecal.graph.nosurgery)
for (i in 1:length(ancom.intersect.nosurgery$taxon)) {
fecal.graph.nosurgery %>%
filter(OTU == ancom.intersect.nosurgery$taxon[i]) %>%
ggplot(aes(x = as.character(Week), y = Abundance, color = Diet)) +
geom_boxplot() + geom_point() +
labs(y="Relative Abundance", title = str_glue(ancom.intersect.nosurgery$taxon[i], " Relative Abundance, fecal"),
x="Week") +
facet_wrap(~Diet) +
theme_bw(base_size = 14) +
theme(legend.position = "top") -> p
print(p)
}
for (i in 1:length(ancom.intersect.nosurgery$taxon)) {
fecal.graph.nosurgery %>%
filter(OTU == ancom.intersect.nosurgery$taxon[i]) %>%
ggplot(aes(x = Diet, y = Abundance, color = Diet)) +
geom_boxplot() + geom_point() +
labs(y="Relative Abundance", title = str_glue(ancom.intersect.nosurgery$taxon[i], " Relative Abundance, fecal"),
x="Diet") +
facet_wrap(~Week) +
theme_bw(base_size = 14) +
theme(legend.position = "top") +
stat_compare_means(method = "wilcox.test", aes(label = after_stat(p.signif)),
label.x = 1.5) -> p
print(p)
}
## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `data.frame()`:
## ! arguments imply differing number of rows: 0, 1
## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `data.frame()`:
## ! arguments imply differing number of rows: 0, 1
DESeq2 provides more information on the results than ANCOM II does in the sense that you can specifically set the contrast. There were a number of taxa that were differentially abundant based on genotype:
fecal.Q2.deseq %>% filter(row %in% tab.fecal.Q2$taxon) %>%
filter(str_detect(Comparison, "Genotype") == TRUE) -> nosurgery.genotype
nosurgery.genotype
fecal.Q2.deseq.intersect.wide %>% filter(row %in% nosurgery.genotype$row)
It looks like these taxa were also affected by diet and week.
for (i in 1:length(nosurgery.genotype$row)) {
fecal.graph.nosurgery %>%
filter(OTU == nosurgery.genotype$row[i]) %>%
ggplot(aes(x = Genotype, y = Abundance, color = Genotype)) +
geom_boxplot() + geom_point() +
labs(y="Relative Abundance", title = str_glue(nosurgery.genotype$row[i], " Relative Abundance, fecal"),
x="Diet") +
facet_grid(vars(Diet), vars(Week)) +
theme_bw(base_size = 14)-> p
print(p)
}
fecal.graph.nosurgery %>%
filter(OTU %in% nosurgery.genotype$row) %>%
group_by(OTU, Diet, Week) %>%
filter(sum(Abundance) > 0) %>% # remove groups that are all zeroes
wilcox_test(Abundance ~ Genotype) %>%
adjust_pvalue(method = "holm") %>%
add_significance()
NS after multiple testing.
fecal.Q1.deseq %>% filter(row %in% ancom.intersect.hfd$taxon) -> fecal.Q1.deseq.intersect
fecal.Q2.deseq %>% filter(row %in% ancom.intersect.nosurgery$taxon) -> fecal.Q2.deseq.intersect
write_csv(fecal.RA.Full, file.path(git.dir, "fecal_relative_abundance.csv"), append = FALSE, col_names = TRUE)
#write_csv(ancom.intersect.hfd, file.path(git.dir, "fecal_ancom_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(ancom.intersect.nosurgery, file.path(git.dir, "fecal_ancom_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)
#write_csv(fecal.Q1.deseq.intersect, file.path(git.dir, "fecal_deseq_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(fecal.Q2.deseq.intersect, file.path(git.dir, "fecal_deseq_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] cowplot_1.1.1 stringr_1.5.0
## [3] tibble_3.2.1 readr_2.1.4
## [5] tidyr_1.3.0 dplyr_1.1.3
## [7] DESeq2_1.40.2 SummarizedExperiment_1.30.2
## [9] Biobase_2.60.0 MatrixGenerics_1.12.3
## [11] matrixStats_1.1.0 GenomicRanges_1.52.1
## [13] GenomeInfoDb_1.36.4 IRanges_2.34.1
## [15] S4Vectors_0.38.2 BiocGenerics_0.46.0
## [17] plotly_4.10.3 DT_0.29
## [19] ggrepel_0.9.4 microbiome_1.22.0
## [21] colorBlindness_0.1.9 RColorBrewer_1.1-3
## [23] pals_1.8 ANCOMBC_2.2.2
## [25] microViz_0.11.0 ggpubr_0.6.0
## [27] ggplot2_3.4.4 rstatix_0.7.2
## [29] fs_1.6.3 phyloseq_1.44.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 DirichletMultinomial_1.42.0
## [3] httr_1.4.7 doParallel_1.0.17
## [5] numDeriv_2016.8-1.1 tools_4.3.1
## [7] doRNG_1.8.6 backports_1.4.1
## [9] utf8_1.2.4 R6_2.5.1
## [11] vegan_2.6-4 lazyeval_0.2.2
## [13] mgcv_1.9-0 rhdf5filters_1.12.1
## [15] permute_0.9-7 withr_2.5.1
## [17] gridExtra_2.3 cli_3.6.1
## [19] sandwich_3.1-0 labeling_0.4.3
## [21] sass_0.4.7 mvtnorm_1.2-3
## [23] proxy_0.4-27 yulab.utils_0.1.0
## [25] foreign_0.8-85 dichromat_2.0-0.1
## [27] scater_1.28.0 decontam_1.20.0
## [29] maps_3.4.1 readxl_1.4.3
## [31] rstudioapi_0.15.0 RSQLite_2.3.3
## [33] generics_0.1.3 gridGraphics_0.5-1
## [35] vroom_1.6.4 gtools_3.9.4
## [37] car_3.1-2 Matrix_1.6-1.1
## [39] biomformat_1.28.0 ggbeeswarm_0.7.2
## [41] fansi_1.0.5 DescTools_0.99.52
## [43] DECIPHER_2.28.0 abind_1.4-5
## [45] lifecycle_1.0.3 multcomp_1.4-25
## [47] yaml_2.3.7 carData_3.0-5
## [49] rhdf5_2.44.0 Rtsne_0.16
## [51] grid_4.3.1 blob_1.2.4
## [53] crayon_1.5.2 lattice_0.22-5
## [55] beachmat_2.16.0 mapproj_1.2.11
## [57] pillar_1.9.0 knitr_1.44
## [59] boot_1.3-28.1 gld_2.6.6
## [61] codetools_0.2-19 glue_1.6.2
## [63] data.table_1.14.8 MultiAssayExperiment_1.26.0
## [65] vctrs_0.6.4 treeio_1.24.3
## [67] Rdpack_2.6 cellranger_1.1.0
## [69] gtable_0.3.4 cachem_1.0.8
## [71] xfun_0.40 rbibutils_2.2.16
## [73] S4Arrays_1.2.1 survival_3.5-7
## [75] SingleCellExperiment_1.22.0 iterators_1.0.14
## [77] gmp_0.7-2 TH.data_1.1-2
## [79] nlme_3.1-163 bit64_4.0.5
## [81] bslib_0.5.1 irlba_2.3.5.1
## [83] vipor_0.4.5 rpart_4.1.21
## [85] colorspace_2.1-0 DBI_1.1.3
## [87] Hmisc_5.1-1 nnet_7.3-19
## [89] ade4_1.7-22 Exact_3.2
## [91] tidyselect_1.2.0 bit_4.0.5
## [93] compiler_4.3.1 htmlTable_2.4.2
## [95] BiocNeighbors_1.18.0 expm_0.999-7
## [97] DelayedArray_0.26.7 checkmate_2.3.0
## [99] scales_1.2.1 digest_0.6.33
## [101] minqa_1.2.6 rmarkdown_2.25
## [103] XVector_0.40.0 htmltools_0.5.6.1
## [105] pkgconfig_2.0.3 base64enc_0.1-3
## [107] lme4_1.1-34 sparseMatrixStats_1.12.2
## [109] fastmap_1.1.1 rlang_1.1.1
## [111] htmlwidgets_1.6.2 DelayedMatrixStats_1.22.6
## [113] farver_2.1.1 jquerylib_0.1.4
## [115] zoo_1.8-12 jsonlite_1.8.7
## [117] energy_1.7-11 BiocParallel_1.34.2
## [119] BiocSingular_1.16.0 RCurl_1.98-1.13
## [121] magrittr_2.0.3 Formula_1.2-5
## [123] scuttle_1.10.3 GenomeInfoDbData_1.2.10
## [125] Rhdf5lib_1.22.1 munsell_0.5.0
## [127] Rcpp_1.0.11 ape_5.7-1
## [129] viridis_0.6.4 CVXR_1.0-11
## [131] stringi_1.7.12 rootSolve_1.8.2.4
## [133] zlibbioc_1.46.0 MASS_7.3-60
## [135] plyr_1.8.9 lmom_3.0
## [137] Biostrings_2.68.1 splines_4.3.1
## [139] multtest_2.56.0 hms_1.1.3
## [141] locfit_1.5-9.8 igraph_1.5.1
## [143] ggsignif_0.6.4 rngtools_1.5.2
## [145] reshape2_1.4.4 ScaledMatrix_1.8.1
## [147] evaluate_0.22 tzdb_0.4.0
## [149] nloptr_2.0.3 foreach_1.5.2
## [151] purrr_1.0.2 rsvd_1.0.5
## [153] broom_1.0.5 Rmpfr_0.9-3
## [155] e1071_1.7-13 tidytree_0.4.5
## [157] viridisLite_0.4.2 class_7.3-22
## [159] gsl_2.1-8 lmerTest_3.1-3
## [161] memoise_2.0.1 beeswarm_0.4.0
## [163] cluster_2.1.4 TreeSummarizedExperiment_2.8.0
## [165] mia_1.8.0